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(3)  Scientific  Progress  and  Accomplishments 


An  Efficient  Electromagnetic-Physics-Based 
Numerical  Technique  For  Modeling  and  Optimization 

of  High-Frequency  Transistors 


Index  Terms—  Adaptive  Grids,  Full  Hydrodynamic  Model,  Global  Modeling,  Maxwell’s 
Equations,  Multifinger  Transistors,  MRTD,  Semiconductor  Simulation,  Wavelets. 

Abstract —  We  present  a  fast  wavelet-based  time-domain  modeling  technique  to  study  the 
effect  of  electromagnetic-wave  propagation  on  the  performance  of  high  power  and 
frequency  multifinger  transistors.  The  proposed  approach  solves  the  active  device  model 
that  combines  the  transport  physics,  and  Maxwell’s  Equations  on  nonuniform  self-adaptive 
grids,  obtained  by  applying  wavelet  transforms  followed  by  hard  thresholding.  This  allows 
forming  line  and  coarse  grids  in  the  locations  where  variable  solutions  change  rapidly  and 
slowly,  respectively.  Comparison  graphs  showed  a  CPU-time  reduction  of  75%  compared 
to  a  uniform-grid  case,  while  maintaining  the  same  degree  of  accuracy.  After  validation, 
the  potential  of  the  developed  technique  is  demonstrated  by  electromagnetic-physical 
modeling  of  multifinger  transistors.  Different  numerical  examples  are  presented; 
emphasizing  that  accurate  modeling  of  high-frequency  devices  should  incorporate  the 
effect  of  EM-wave  propagation  and  electron-wave  interactions,  within  and  around  the 
device.  Moreover,  high-frequency  advantages  of  multifinger  transistors  over  single-finger 
transistors  are  highlighted.  To  our  knowledge,  this  is  the  first  time  in  literature  to 
introduce  and  implement  a  fully  numerical  electromagnetic-physics-based  simulator  for 
accurate  modeling  of  high-frequency  transistors. 


I.  Introduction 


MUETIFINGER  transistors  have  proven  better  performance  over  conventional  transistors, 
especially  at  very  high  frequency  [1]-[11].  However,  till  now,  modeling  of  such  devices 
did  not  account  for  EM-wave  effects  as  well  as  electron-wave  interactions  using  a  fully 
numerical  simulator.  Accordingly,  it  is  indispensable  to  present  analysis  of  multi  linger 
transistors  based  on  a  coupled  electromagnetic-physics-based  simulator.  The  possibility  of 
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achieving  this  type  of  modeling  is  addressed  by  global  circuit  modeling  that  has  been 
demonstrated  in  [12] -[17]. 

Global  modeling  is  a  tremendous  task  that  involves  advanced  numerical  techniques  and 
different  algorithms.  As  a  result,  it  is  computationally  expensive  [17].  Therefore,  there  is  an 
imperative  need  to  develop  a  new  approach  to  reduce  the  simulation  time,  while  maintaining  the 
same  degree  of  accuracy  presented  by  global  modeling  techniques.  One  approach  is  to 
adaptively  refine  grids  in  locations  where  the  unknown  variables  vary  rapidly.  Such  technique  is 
called  multiresolution  time  domain  (MRTD),  and  a  very  attractive  way  to  implement  it  is  to  use 
wavelets  [18]-[19]. 

MRTD  approach  has  been  successfully  applied  to  finite-difference  time-domain  (FDTD) 
simulations  of  passive  structures  [20]-[31].  However,  for  the  active  devices,  that  are 
characterized  by  a  set  of  coupled  and  highly  nonlinear  partial  differential  equations,  applying  the 
same  approach  would  become  quite  time  consuming  [32],  Several  different  approaches  for 
solving  partial  differential  equations  (PDE’s)  using  wavelets  have  been  considered.  It  has  been 
observed  by  several  authors  that  nonlinear  operators  such  as  multiplication  are  too 
computationally  expensive  when  conducted  directly  on  a  wavelet  basis.  One  of  the  approaches 
for  solving  PDE’s  is  the  Interpolating  Wavelets  technique  presented  in  [33],  in  which  the 
nonlinearities  are  dealt  with  using  the  so-called  sparse  point  representation  (SPR).  Interpolating 
Wavelets  have  been  successfully  applied  to  the  simple  drift  diffusion  active  device  model  [34]- 
[36].  Being  primarily  developed  for  long-gate  devices,  the  drift  diffusion  model  leads  to 
inaccurate  estimations  of  device  internal  distributions  and  microwave  characteristics  for 
submicrometer  devices  [37].  It  is  worth  mentioning  that  in  [33],  the  author  proposed  a  new 
technique  to  solve  simple  forms  of  Hyperbolic  PDE’s  using  an  Interpolating  Wavelet  scheme. 
These  PDE’s  can  represent  Maxwell’s  Equations  or  the  simple  drift-diffusion  model  but  not  the 
complete  hydrodynamic  model.  Thus,  a  new  approach  to  apply  wavelets  to  the  hydrodynamic 
model  PDE’s  is  needed,  along  with  extending  it  to  Maxwell’s  Equations,  for  accurate  modeling 
of  submicrometer  devices,  while  achieving  a  CPU -time  reduction. 

In  this  paper,  a  unified  approach  to  apply  wavelets  to  the  full  hydrodynamic  model  and 
Maxwell’s  equations  is  developed.  The  main  idea  is  to  take  snapshots  of  the  solution  during  the 
simulation,  and  apply  wavelet  transfonn  to  the  current  solution  to  obtain  the  coefficients  of  the 
details.  The  coefficients  of  the  details  are  then  normalized,  and  a  threshold  is  applied  to  obtain  a 
nonuniform  grid.  Two  independent  grid-updating  criteria  are  developed  for  the  active  and 
passive  parts  of  the  problem.  Moreover,  a  threshold  formula  that  is  dependent  on  the  variable 
solution  at  any  given  time  has  been  developed  and  verified.  In  addition,  a  full-wave  global 
modeling  simulator  is  developed  to  study  the  EM-wave  propagation  effect  on  high  power  and 
frequency  multifinger  transistors.  A  comprehensive  set  of  results  is  included  along  with 
illustrative  comparison  graphs. 

This  paper  is  organized  as  follows.  Section  II  reviews  briefly  the  theory  of  MRTD.  Problem 
description  is  given  in  Section  III.  Full  descriptions  of  the  proposed  algorithm  along  with 
illustrative  graphs  are  provided  in  Section  IV.  Technique  validation  is  presented  in  Section  V. 
Error  and  stability  analysis  are  discussed  in  Section  VI.  While,  the  microwave  characteristics  of 
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high-frequency  transistors  are  provided  in  Section  VII,  and  results  of  electromagnetic  physical 
modeling  of  multifinger  transistors  are  presented  in  Section  VIII.  Finally,  conclusions  are 
provided  in  Section  X. 


II.  Fundamentals  Of  MRTD 

The  construction  of  biorthogonal  wavelet  bases  relies  on  the  notation  of  multiresolution 
analysis  [38].  This  notation  gives  a  formal  description  of  the  intuitive  idea  that  every  signal  can 
be  constructed  by  a  successive  refinement,  by  iteratively  adding  details  to  an  approximation.  The 
coefficients  of  the  approximations  are  given  by: 

+oo 

ax[n,m\=  I  x{t)<pm{t)dt  (1) 

—00 

where  cpnm  ( t )  is  the  family  of  dilates  and  translates  of  the  scaling  function  formed  as: 

*>„(<)  =  2"'>(2-<-«).  U) 

On  the  other  hand,  Eq.  3  gives  the  coefficients  of  the  details. 

+00 

dx[n,m]=  j  x(t)y/nm(t)dt.  (3) 

—  00 

Where  y/nm  (t)  is  the  family  of  dilates  and  translates  of  the  wavelet  function  defined  as: 

Wnmit)  =  2m/2¥(2mt-n);n,m^Z}.  (4) 

While  some  wavelets  such  as  Daubechies  are  asymmetrical  [38],  it  is  possible  to  create 
symmetric  wavelets  with  compact  support  by  using  two  sets  of  wavelets,  one  to  compose  the 
signal  and  the  other  to  construct  it.  Such  wavelets  are  called  Biorthogonal  [39]. 

III.  Problem  Description 

The  transistor  model  used  in  this  work  is  a  3D  full-hydrodynamic  large-signal 
electromagnetic-physical  model.  The  active  device  model  is  based  on  the  moments  of  the 
Boltzmann’s  Transport  Equation  obtained  by  integrating  over  the  momentum  space.  The 
integration  results  in  a  strongly  coupled  highly  nonlinear  set  of  partial  differential  equations, 
called  the  conservation  equations.  These  equations  provide  a  time-dependent  self-consistent 
solution  for  carrier  density,  carrier  energy,  and  carrier  momentum,  which  are  given  as  follows. 

■  current  continuity 

■  energy  conservation 


dn 

dt 


+  V.(nu)  =  0. 


d(ne) 

8t 


+  qnv.E  +  V.(w(s  +  KBT)  )  =  - 


fl(g-So) 

tJe) 


(5) 

(6) 
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x-momentum  conservation 


niPx-Po) 

*M00 


(7) 


d(npx) 

dt 


+  clnEx  +xr(nP,ux  +  nKBT )  = 

ox 


In  the  above  equations,  n  is  the  electron  concentration,  o  is  the  electron  velocity,  E  is  the 
electric  field,  s  is  the  electron  energy,  s0  is  the  equilibrium  thennal  energy,  and  p  is  the  electron 

momentum.  The  energy  and  momentum  relaxation  times  are  given  by  r£  and  xm  (  respectively. 

Similar  expression  is  obtained  for  the  y-direction  momentum.  The  three  conservation  equations 
are  solved  in  conjunction  with  Maxwell’s  Equations: 


VxE 

V  x  H 


8B 

dt 


(I) 

dt 


+  J 


(8) 

(9) 


where  E  is  the  electric  field,  H  is  the  magnetic  field,  D  is  the  electric  flux  density,  and  B  is 
the  magnetic  flux  density.  The  fields  in  Maxwell’s  Equations  are  updated  using  the  current 
density  /estimated  by  Eq.  (10). 


J(t)  =  —qnu(t).  (10) 

The  low  field  mobility  is  given  by  the  empirical  relation  [40]. 

,  8000  17,..icmi/V.s.  (11) 

The  above  model  accurately  describes  all  the  non- stationary  transport  effects  by 
incorporating  energy  dependence  into  all  the  transport  parameters  such  as  effective  mass  and 
relaxation  times,  along  with  including  EM-wave  effects.  Fig.  1  shows  the  cross-section  of  the 
simulated  structure  used  to  validate  the  new  algorithm,  with  parameters  summarized  in  Table  I. 


I 

y 


Fig 


0.5  |Jm 

+ - ► 

Source 


03  [im 
Gate 


Drain 


Doped  GaAs  (2x10' 7  cnO 


03  [im 


Undoped  GaAs 


2B  [im 


1.  Cross-section  of  the  simulated  transistor. 
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TABLE  I 


Transistor  Parameters  Used  In  The  Simulation 


Drain  and  source  contacts 

0.5  pm 

Gate-source  separation 

0.5  pm 

Gate-drain  separation 

1.0  pm 

Device  thickness 

0.8  pm 

Device  length 

2.8  pm 

Gate  length 

0.3  pm 

Device  Width 

250  pm 

Active  layer  thickness 

0.2  pm 

Active  layer  doping 

2xl017  cm"3 

Schottky  barrier  height 

0.8  V 

DC  gate-source  voltage 

-0.5  V 

DC  drain-source  voltage 

3.0  V 

IV.  The  Proposed  Algorithm 


Fig.  2  shows  the  flow  chart  of  the  proposed  algorithm.  A  uniform  grid  is  defined  at  the 
beginning  of  the  simulation.  Equations  5  through  7  are  then  solved  in  the  sequence  shown  by  the 
flow  chart  to  update  the  grid  of  the  different  variables  at  the  new  iteration  with  the  following 
criterion: 


x 


/ 

m  ax  ,m  in 


-  X 


c 

m  ax  ,m  in 


X 


l 

m  ax  ,m  in 


>  8 


H  y  d  ro 


(12) 


The  updating  criterion  checks  if  the  solution  of  the  variable  x  has  changed  by  SHydro  since  last 

iteration  using  wavelet  transform.  The  subscripts  c  and  /  designate  quantities  defined  in  the 
current  time  and  last  time  where  wavelet  transform  is  performed,  respectively.  The  subscript 
“max,  min”  indicates  that  the  maximum  and  minimum  of  the  variable  x  are  checked  with 
Relation  12  at  the  same  time.  It  is  significant  to  note  here  that  boundary  grid-points  are  not 
included  for  the  maximum  or  minimum  checking.  The  value  of  8Uydro  used  in  the  simulation  is 

0.1.  If  Relation  12  is  satisfied,  wavelet  transform  is  performed  on  the  current  variable  solution 
followed  by  hard  thresholding  to  obtain  a  new  nonuniform  grid  for  the  variable  x.  Biorthogonal 
wavelets  are  used  with  notation  BI03.1  to  point  out  three  vanishing  moments  for  the  mother 
wavelet  and  only  one  vanishing  moment  for  the  scaling  function.  The  nonuniform  grids  of  the 
different  variables  are  then  combined  into  only  one  nonuniform  grid  for  the  next  iteration.  The 
above  steps  are  repeated  until  the  stopping  criterion  is  satisfied. 
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Fig.  2.  Generic  flow  chart  of  the  proposed  algorithm. 
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It  should  be  noted  that  magnitude  ranges  of  the  variables  used  in  the  simulations  vary 
dramatically.  For  instance,  carrier  density  per  cm  ,  is  on  the  order  of  10  ,  while,  energy 
expressed  in  eV ,  is  on  the  order  of  0.5 .  Accordingly,  the  threshold  value  should  be  dependent 
on  the  variable  solution  at  any  given  iteration.  The  proposed  threshold  formula  is  given  as: 


In  this  equation,  T0  is  the  initial  threshold  value,  ci  ’  s  are  the  coefficients  of  the  details,  and 

N  is  the  number  of  grid  points  in  the  x-  or  y-direction.  Hence,  the  value  of  the  threshold  T 
depends  mainly  on  the  variable  solution  at  any  given  time  rather  than  being  fixed.  The  values  of 
T0  used  in  the  simulation  are  0.001,  0.01,  and  0.05,  respectively. 

It  is  significant  to  note  that  the  proposed  algorithm  is  carefully  developed  such  that  it  is 
general  and  device  independent.  This  includes  the  threshold  value  given  by  Eq.  13.  Considering 
this  equation,  it  is  clear  that  it  has  no  device-dependent  parameters.  Furthermore,  the  different 
values  of  T0  are  used  only  to  investigate  the  tradeoffs  between  CPU-time  and  accuracy,  which  is 
a  general  approach  provided  by  papers  introducing  similar  threshold-controlled  wavelet-based 
techniques,  for  instance  in  [21]  and  [35], 

In  this  paper,  a  new  technique  to  conceive  the  nonunifonn  grids  using  wavelets  has  been 
developed.  The  main  idea  is  to  apply  wavelet  transform  to  the  variable  solution  at  any  given 
time  to  obtain  the  coefficients  of  the  details,  which  are  then  normalized  to  its  maximum.  Only 
grid  points  where  the  value  of  the  normalized  coefficients  of  the  details  larger  than  the  threshold 
value  given  by  Eq.  13  are  included.  Fig.  3  exemplifies  how  the  proposed  algorithm  obtains  the 
nonuniform  grid  using  longitudinal  compression  only.  For  instance,  Fig.  3. a.  shows  the 
nonnalized  amplitude  of  the  coefficients  of  the  details  for  the  electron  energy.  While,  Fig.  3.b. 
marks  the  grid  points  remaining  after  thresholding  the  normalized  coefficients  of  the  details 
using  Eq.  13.  It  should  be  observed  that  the  proposed  technique  accurately  removes  grid  points 
in  the  locations  where  variable  solutions  change  very  slowly. 

The  overall  grid  obtained  needs  further  processing  in  order  to  define  a  finite-difference 
scheme  on  it.  The  simplest  way  to  achieve  that  is  to  have  the  same  number  of  grid  points  for  the 
parallel  cross-sections,  while  the  number  of  grid  points  in  the  longitudinal  cross-sections  and  the 
transverse  cross-sections  need  not  to  be  the  same.  Following  the  above  procedure,  it  was  found 
that  boundary  conditions  implementation,  including  Ohmic  and  Schottky  contacts,  does  not  need 
special  treatment.  They  can  be  treated  similar  to  the  standard  finite-difference  (FD)  scheme.  The 
only  issue  the  algorithm  needs  to  keep  track  of  is  identifying  the  new  boundaries  of  the  metallic 
contacts  for  each  new  grid,  which  is  straightforward. 
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Electron  Energy  (ev)  Normalized  Amplitude 


(b)  (d) 

Fig.  3.  (a)  Normalized  details  coefficients  for  the  electron  energy  at  a  certain  longitudinal  cross-section,  (b) 
Grid  points  marked  on  the  actual  curve  for  the  electron  energy  at  the  same  longitudinal  cross-section,  (c) 
Normalized  details  coefficients  for  the  x-momentum  at  a  certain  longitudinal  cross-section,  (d)  Grid  points 
marked  on  the  actual  curve  for  the  x-momentum  at  the  same  longitudinal  cross-section. 
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TABLE  II 


Grid  Adaptability  of  the  Different  Variables  For  T()  =  1% 


Unknowns  Remaining  After  Unknowns  Remaining  After 

Total 

Variable 

Transverse  Compression  Longitudinal  Compression 

Unknowns 

(%) 

Time  Iteration#  120 

(%) 

Remaining 

(%) 

Potential 

5.69 

7.74 

0.63 

Carrier  Density 

6.54 

14.92 

2.64 

Energy 

39.65 

17.63 

8.54 

x  -Momentum 

43.39 

16.06 

9.18 

y  -Momentum 

16.11 

17.53 

3.78 

All  Variables 

65.14 

Time  Iteration  #  250 

22.36 

14.43 

Potential 

5.88 

8.59 

0.76 

Carrier  Density 

13.69 

16.70 

5.18 

Energy 

39.21 

23.00 

12.26 

x  -Momentum 

43.65 

19.09 

10.64 

y  -Momentum 

20.02 

19.46 

7.95 

All  Variables 

61.94 

Time  Iteration  #  480 

28.93 

20.58 

Potential 

6.27 

9.23 

0.90 

Carrier  Density 

21.51 

17.16 

7.71 

Energy 

43.99 

28.88 

15.72 

x  -Momentum 

38.57 

23.44 

12.72 

y  -Momentum 

26.20 

26.76 

13.53 

All  Variables 

58.84 

Time  Iteration  #  590 

36.25 

25.05 

Potential 

6.04 

9.64 

0.93 

Carrier  Density 

29.88 

18.24 

10.61 

Energy 

48.85 

31.88 

18.43 

x  -Momentum 

41.91 

29.08 

16.21 

y  -Momentum 

32.91 

37.04 

20.36 

All  Variables 

62.36 

Time  Iteration  #  730 

44.73 

31.74 

Potential 

7.01 

11.13 

1.15 

Carrier  Density 

34.08 

16.55 

8.42 

Energy 

39.77 

35.64 

18.43 

x  -Momentum 

41.91 

27.98 

13.32 

y  -Momentum 

51.46 

34.84 

25.00 

All  Variables 

62.84 

58.96 

36.43 

Table  II  shows  the  evolution  of  the  nonuniform  grids.  It  can  be  observed  that  the  number  of 
grid  points  for  the  overall  grid  increases  as  time  advances.  The  reason  is  that  at  the  beginning  of 
the  simulation  the  solution  is  not  completely  fonned  yet,  and  as  time  marches,  more  grid  points 
are  needed  to  incorporate  the  changes  in  the  solution.  Furthermore,  the  different  variable  grids 
should  not  be  updated  at  the  same  rate.  For  instance,  it  is  apparent  that  the  potential  needs  not  to 
be  updated  at  the  same  rate  as  the  other  variables.  Notice  that  Table  II  is  used  for  illustration 
purposes  to  demonstrate  the  way  the  different  variable  grids  change.  In  the  actual  simulator,  the 
potential  grid  is  updated  a  few  times  at  the  beginning  of  the  simulation,  and  then  it  remains 
unchanged. 
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Now,  we  turn  our  attention  to  Maxwell’s  Equations.  The  passive  part  of  the  FET  represents  a 
co-planar  structure  in  which  a  3D  FDTD  is  developed  to  solve  for  the  electric  and  magnetic 
fields.  The  current  density  estimated  from  the  active  device  conversation  equations  is  used  to 
update  the  variables  in  Maxwell’s  Equations. 

It  is  importance  to  state  that  the  same  approach  developed  to  obtain  the  nonuniform  grid  for 
the  variables  of  the  conversation  equations  is  applied  to  Maxwell’s  Equations  as  well.  However, 
a  different  updating  mechanism  should  be  developed  to  keep  track  of  the  wave  propagation 
within  the  passive  structure.  The  following  is  the  algorithm  developed  for  the  grid  updating  of 
FDTD  simulations. 


Step  1 :  Construct  a  3D  matrix  M  that  has  only  0’s  and  1  ’s,  based  whether  or  not  we  have  a  non¬ 
zero  solution  of  the  field  at  this  location.  For  example,  “1”  is  assigned  if  a  non-zero  field  solution 
exists,  and  “0”  elsewhere. 

Step  2:  Estimate  the  value  of  8jim)  (FDTD  grid-updating  factor)  as: 


8 


FDTD 


y  (m  ©  m  , , ) . . , 

v  new  old  7  i,j,k 

i,j,k _ 

N xd^ ydN _d 


(14) 


Where  Mnew  and  Mold  are  the  matrices  constructed  using  step  one  for  the  current,  and  old 
solutions  of  the  fields,  respectively.  Nxd ,  Nyd  ,  and  Nzd  are  the  number  of  grid  points  in  x,  y  , 
and  z  directions,  respectively. 

Step  3:  Check  8m[D,s  value  against  a  predefined  value,  for  example  5%. 

Step  4:  If  satisfied,  move  the  grid  to  z  =  z  +  dz  .  Where  dz  is  proportional  to  ^fdtd- 
Step  5:  t  =  t  +  dt 


Figure  4  illustrates  examples  of  how  the  nonuniform  grids  are  obtained  for  the  magnetic  and 
electric  fields  at  a  specific  cross-section  for  FDTD  simulations.  For  instance,  Fig.  4. a. shows  the 
nonnalized  amplitude  of  the  coefficients  of  the  details  for  the  electric  field.  While  Fig.  4.b. 
marks  the  grid  points  remaining  after  thresholding  the  normalized  coefficients  of  the  details 
using  Eq.  13.  It  is  observed  that  the  proposed  technique  accurately  removes  the  grid  points  in 
locations  where  variable  solutions  change  very  slowly.  This  would  have  an  effect  of  reducing 
the  CPU-time  by  removing  the  redundant  grid-points  introduced  by  the  original  fonnulation.  The 
total  grid  of  the  electric  field  can  be  achieved  by  obtaining  two  separate  grids  for  the  transverse 
and  longitudinal  compressions,  respectively.  Then  the  two  grids  are  combined  together  using 
logical  ‘AND’  to  conceive  the  overall  grid  for  the  electric  field  at  this  given  time.  It  is  worth 
mentioning  here  that  the  excitation  wave  exists  at  the  source  plan  at  all  times,  and  the  technique 
proposed  here  is  generic  that  can  be  applied  also  to  a  short  pulse  propagating  in  the 
computational  domain. 
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(a) 


(b) 


Fig.  4.  Demonstration  of  the  procedure  employed 
and  magnetic  fields  for  FDTD  simulations. 


magnetic  field  (Amp. /cm) 


obtain  the  nonuniform  grid  for  the  y-direction  electric 


V.  Technique  Validation 
A.  Hydrodynamic  Model  Simulation  Results 

The  approach  presented  in  this  paper  is  general  and  can  be  applied  to  any  unipolar  transistor. 
To  demonstrate  the  potential  of  this  approach,  it  is  applied  to  an  idealized  FET  structure,  which 
is  discretized  by  a  mesh  of  64  Ar  by  64  Av  with  At  =  0.001  ps .  Forward  Euler  is  adopted  as  an 
explicit  finite-difference  method.  In  addition,  upwinding  is  employed  to  have  a  stable  finite- 
difference  scheme.  The  space  step  sizes  are  adjusted  to  satisfy  Debye  length,  while  the  time  step 
value  A t  is  chosen  to  satisfy  the  Courant-Friedrichs-Levy  (CFL)  condition.  First,  DC  simulations 
are  performed  following  the  flow  chart  given  by  Fig.  2,  and  the  current  density  is  calculated 
using  Eq.  10.  DC  excitation  is  performed  by  forcing  the  potential  to  be  equal  to  the  applied 
voltages  to  the  electrodes  (i.e.,  Dirichlet  boundary  conditions). 

It  is  important  to  note  that  a  suitable  approach  to  investigate  the  capabilities  of  the  proposed 
technique  is  to  compare  it  to  the  uniform-grid  algorithm.  In  this  case,  the  new  simulator  will  be 
accurately  evaluated.  Since  both  algorithms,  the  wavelet-based  and  the  uniform  one,  will  run  on 
the  same  computer.  In  addition,  both  algorithms  will  have  the  same  discritization  schemes  and 
the  exact  semiconductor  parameters. 

Fig.  5  shows  the  drain  current  convergence  curves  versus  the  CPU-time  in  seconds  for  the 
cases  of  the  uniform  grid  and  the  proposed  wavelet-based  adaptive  grids  with  different  initial 
threshold  values  T0.  Fig.  5  demonstrates  that  using  the  proposed  wavelet-based  grids  approach 

reduces  the  CPU-time  dramatically.  For  instance,  there  is  a  reduction  of  about  75%  in  CPU-time 
over  the  uniform  grid  case  for  the  initial  threshold  value  of  1%,  while  the  DC  drain  current  error 
is  within  1%.  In  addition,  increasing  the  initial  threshold  beyond  certain  value  has  a  negative 
effect  on  the  accuracy  of  the  final  results.  This  is  apparent  for  T0  equals  to  5%,  where  there  is  no 
agreement  between  the  results  achieved  using  the  uniform  grid  case  and  the  wavelet-based 
nonuniform  grids.  The  reason  is  that  using  large  values  of  T0  implies  that  more  grid  points  are 

removed,  including  important  grid  points  that  will  have  a  negative  effect  on  the  final  result.  On 
the  other  hand,  using  a  very  small  threshold  values  implies  redundant  grid  points.  In  summary, 
there  should  be  an  optimal  value  of  T0  such  that  both  the  CPU-time  and  error  are  minimized.  In 

this  work,  T0  of  1%  is  suggested  to  have  a  considerable  reduction  in  CPU-time,  while  keeping 
error  within  an  acceptable  range. 

The  problem  presented  in  this  paper  is  simply  to  validate  the  new  algorithm,  and  to 
emphasize  that  practical  problems  have  solutions  that  change  very  rapidly  only  in  specific 
domains.  Accordingly,  an  algorithm  could  be  implemented  to  exploit  that  by  solving  the 
equations  on  multiresolution-nonuniform  grids,  obtained  using  wavelets.  Obviously,  the  amount 
of  CPU-time  reduction  depends  on  the  problem  under  consideration.  More  or  less  reduction 
could  be  achieved  for  other  types  of  problems  depending  on  the  sharpness  and  distribution  of 
solutions  in  the  computational  domain. 
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Fig.  5.  DC  drain  current  convergence  curves  for  the  uniform  grid  and  the 
proposed  wavelet-based  non-uniform  grids  for  different  initial  threshold 
values. 


B.  FDTD  Simulation  Results 

A  3D  Yee-based  FDTD  code  is  developed,  with  the  proposed  algorithm  employed.  A 
Guassian  excitation  pulse  is  applied  to  evaluate  the  algorithm  over  a  wide  range  of  frequencies. 
Table  II  shows  that  as  threshold  value  increases,  CPU-time  and  error  introduced  decreases  as 
well.  It  is  noteworthy  to  point  out  that  using  an  initial  threshold  value  equals  to  10%  seems  to 
reduce  error  along  with  the  CPU-time.  However,  considering  Fig.  6,  one  should  conclude  that 
using  T0  equals  to  10%  introduces  dispersion,  which  is  a  serious  type  of  error.  Accordingly,  an 

initial  threshold  value  of  5%  is  recommended  in  terms  of  both  CPU-time  and  error  for  FDTD 
simulations.  It  is  important  to  emphasize  that  the  passive  and  active  parts  of  the  problem  have 
different  optimal  threshold  values.  This  is  expected  since  the  variables  in  the  conservations 
equations  are  highly  nonlinear  compared  to  the  fields  obtained  when  solving  Maxwell’s 
Equations.  Research  leading  to  the  work  presented  in  this  paper  can  be  found  in  [41] -[43], 
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TABLE  III 

Effect  of  the  Threshold  Value  on  Error  and  CPU-Time  For  FDTD  Simulations 


T 

0.0  (Uniform  Grid) 

CPU-Time 

(Seconds) 

744.90 

Error 

2-norm 

co- norm 

0.1% 

300.17 

0.0873% 

8.80% 

1.0% 

205.92 

0.0871% 

8.75% 

5.0% 

155.10 

0.0778% 

7.69% 

10.0% 

111.05 

0.0473% 

3.66% 

Fig.  6.  Potential  of  the  gate  at  a  specific  cross-section  versus  time  for  the  uniform 
grid  case  and  the  proposed  MRTD  algorithm  with  different  values  of  T0 . 
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VI.  Scheme  Error  and  Stability  Analysis 


It  is  important  to  mention  here  that  the  simulation  and  physical  times  are  completely  separate 
entities.  The  simulation  time  required  to  model  a  specific  physical  process  should  vary 
depending  on  the  technique  implemented  in  the  simulation. 

The  purpose  of  this  section  is  to  demonstrate  that  the  mechanism  by  which  error  is 
introduced  when  employing  the  proposed  wavelet-based  technique  is  different  than  of  the 
uniform-grid  case.  The  local  truncation  error  for  the  uniform  grid  case  is  dependent,  in  general, 
on  the  mesh  spacing  ( Ar  and  At7 ),  and  the  time  step  used  At .  On  the  other  hand,  the  local 
truncation  error  for  the  wavelet-based  nonuniform  grids  approach  depends  on  how  accurately  the 
important  grid  points  are  reserved  as  well  as  the  time  step  used.  This  suggests  that  the  local 
truncation  errors,  due  to  spatial  discritization,  for  the  uniform  grid  case  and  the  wavelet-based 
nonuniform  grids  are  different.  The  local  truncation  error  accumulates  from  iteration  to  iteration. 
The  total  truncation  or  discritization  error  is  thus  dependent  on  the  number  of  iterations  used 
(space  and  time  iterations  combined).  Accordingly,  one  can  conclude  that  the  total  error 
introduced  by  the  wavelet-based  technique  due  to  the  local  discritization  errors  accumulating 
during  the  simulation,  may  or  may  not  be  larger  than  that  of  the  uniform  grid  case,  at  least  for  the 
two  cases  of  T0  =0.1%  and  T0=  1%  .  The  reason  is  the  number  of  iterations  required  reaching  the 

steady  state  solution  for  the  uniform  grid  case  is  much  larger  than  that  of  the  proposed  algorithm. 
In  summary,  the  total  error  introduced  depends  on  the  local  truncation  error  along  with  the 
number  of  iterations  required  to  reach  the  final  solution.  This  explains  the  results  in  the  paper 
comparison  figures,  where  it  would  be  difficult  to  draw  a  precise  conclusion  of  which  technique 
is  more  accurate.  This  is  because  for  each  case  or  curve,  the  number  of  iterations  required  to 
obtain  the  steady  state  solution  and  the  local  discritization  errors  are  different.  The  problem  of 
identifying  the  most  accurate  solution  becomes  even  more  difficult  since  we  are  dealing  with  a 
highly  nonlinear  problem. 

It  is  significant  to  call  attention  to  the  fact  that  the  proposed  algorithm  does  not  have  any 
stability  constraints  if  At  is  chosen  to  satisfy  the  CFL  condition  at  the  beginning  of  the 
simulation.  The  reason  is,  as  the  simulation  progresses,  the  spatial  distances  employed  become 
even  larger  than  the  ones  introduced  at  the  beginning.  This  represents  an  extra  benefit  of  using 
the  proposed  algorithm  that  it  does  not  need  any  time-step  At  change  while  the  simulation  is  in 
progress. 
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VII.  Microwave  Characteristics  of  High-Frequency  transistors 


To  study  the  characteristics  of  transistors  at  high  frequency,  a  time-domain  Guassian  signal 
is  applied  between  the  source  and  gate  electrodes.  The  input  and  output  time-domain  signals  are 
observed  at  different  points  along  the  width  of  the  device.  The  characteristics  of  the  device  are 
then  estimated.  For  example,  the  propagation  constant  y  can  be  evaluated  as: 


y  =  I|oge 


\F(co,z+l)\ 
1  F(cd,z)  /• 


(15) 


Where  F(co,z )  is  the  Fourier  transform  of  the  time-domain  signal.  The  attenuation  and 

propagation  constants  are  evaluated  as  the  real  and  imaginary  parts  of  y,  respectively.  Fig.  7 
shows  the  attenuation  constant  as  a  function  of  frequency  at  different  points  along  the  device 
width.  Considering  Fig.  7,  it  should  be  noticed  that  the  attenuation  constant  increases  with 
frequency  as  well  as  from  point  to  point  along  the  device  width. 


Fig.  7.  Attenuation  constant  as  a  function  of  frequency  at  different  points 
along  the  device  width  for  the  gate  electrode. 
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The  phase  velocity  v  h  and  effective  dielectric  constant  s;.  can  be  estimated  using  Equations 
(16)-(17),  respectively. 


v 


ph 


CO_ 

P 


(16) 


P2c 

( 2/r/)“ 


Where  P  is  the  propagation  constant,  c  is  the  free-space  wave  velocity,  and  co  is  the 
frequency  in  rad./sec.  Figures  (8)-(9)  show  the  effective  dielectric  constant  and  phase  velocity 
versus  frequency  at  different  points  along  the  device  width,  respectively.  The  results  shown  in 
Figures  (8)-(9)  are  mainly  due  to  the  change  of  the  distribution  of  the  electric  field  as  a  function 
of  frequency  and  distance.  The  results  presented  in  Figures  (8)-(9)  coincide,  conceptually,  with 
those  presented  in  [44].  These  results  are  distinctive  to  high  frequency  devices  only,  which  could 
be  minimized  by  employing  optimized  microwave  structures  such  as  multifinger  transistors. 


Fig.  8.  Effective  dielectric  constant  as  a  function  of  frequency  at  different 
points  along  the  device  width  for  the  gate  electrode. 
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Fig.  9.  Phase  velocity  as  a  function  of  frequency  at  different  points  along 
the  device  width  for  the  gate  mode. 

VIII.  EM-wave  Effects  on  High-Frequency  Multifinger  Transistors 

A.  EM-  Wave  Propagation  effects  and  Electron-  Wave  Interactions 

In  this  section,  a  full-wave  physical  simulator  is  developed  to  model  two  closely  packed 
millimeter- wave  transistors.  Fig.  10  gives  a  3D  view  of  the  simulated  transistors.  The  simulated 
devices  are  biased  to  V/s  =  3.0V  and  Vos  =  -0.1V .  The  gate-length  for  the  transistors  is  set  to 

0.2 pm  .  The  DC  distributions  are  obtained  by  solving  the  active  device  model  only.  A  sinusoidal 
signal  is  employed  in  the  AC  simulations  with  peak  value  of  lOOmV  and  frequency  of  80  GHz, 
respectively.  The  two  transistors  shown  in  Fig.  10  are  identical.  First,  full- wave  simulations  are 
carried  out  for  one  transistor  only,  and  the  results  are  depicted  in  Fig.  1 1 .  Considering  this  figure, 
one  should  observe  the  variations  of  the  output  voltage  with  distance  along  the  device  width.  The 
reasons  are  due  to  the  nonlinear  energy  build-up  along  the  device  width,  and  due  to  the  phase 
velocity  mismatch  between  the  EM-waves  at  the  gate  and  drain  electrodes.  Fig.  1 1  demonstrates 
the  importance  of  coupling  the  EM-waves  with  the  semiconductor  transport  physics  for  accurate 
modeling  of  millimeter- wave  transistors. 

Now,  we  turn  our  attention  to  the  full-wave  simulation  of  the  two  transistors  shown  in  Fig. 
10.  First,  we  assume  that  one  of  the  transistors  is  operating,  while  the  other  transistor  is  not.  Fig. 
12  depicts  the  simulation  results,  which  emphasize  the  significance  to  include  the  EM-wave 
propagation  effects,  not  only  inside  the  device,  but  around  it  as  well.  In  fact,  this  is  the  basic 
theory  of  operation  of  multifinger  transistors.  Ideally,  the  non-operating  transistor  should  have  a 
zero  drain  potential,  however  due  to  the  proximity  of  an  operating  transistor,  an  induced  voltage 
that  varies  along  the  device  width  is  introduced. 
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Fig.  10.  3D  view  of  the  simulated  transistors  ( not  to  scale). 


Next,  the  two  transistors  in  the  configuration  shown  in  Fig.  10  are  simulated,  assuming 
that  both  transistors  are  now  operating.  There  may  be  two  cases  to  consider.  The  first  case  is  to 
assume  the  drains  of  the  two  transistors  are  adjacent  to  each  other  (the  case  of  multifinger 
transistors).  While,  the  other  case  is  to  consider  the  drain  of  one  of  the  transistors  is  adjacent  to 
the  source  of  the  other  transistor.  Figures  (13)-(14)  show  the  simulation  results.  The  first 
conclusion  that  can  be  drawn  out  of  the  two  figures  is  that  the  proximity  of  an  operating 
transistor  affects  the  output  voltage  due  to  the  EM-wave  propagation.  Furthermore,  the  EM-wave 
effects  for  the  case  of  two  adjacent  drain  electrodes  is  much  larger  than  the  other  case.  This  is 
expected,  since  the  drain  electrode  has  the  amplified  output  signal.  It  is  important  to  mention  that 
the  results  in  Figures  (11),  (13),  and  (14)  are  normalized,  which  do  not  represent  the  actual  gain. 


Time  (s)  x  10  " 

Fig.  11.  Drain  voltage  (normalized)  of  the  simulated  transistor  when 
EM-wave  propagation  and  electron-wave  interaction  are  considered  at 
different  points  in  the  z-direction. 


24 


Fig.  12.  Potential  of  a  passive  electrode  at  different  points  in  the  z-direction 
induced  due  to  the  proximity  of  an  operating  transistor  excited  by  a  Gaussian 
signal. 


Fig.  13.  Drain  voltage  (normalized)  at  z  =  62.5///??  when  EM-wave 
propagation  and  electron-wave  interaction  are  considered.  Solid  line: 
transistor  is  simulated  alone.  Dashed  line:  Source  electrode  of  a  second 
operating  transistor  is  0.5 /um  apart  from  the  drain  of  the  simulated  transistor. 
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Fig.  14.  Drain  voltage  (normalized)  at  z  =  62.5//7J7  when  EM-wave 
propagation  and  electron-wave  interactions  are  considered.  Solid  line: 
transistor  is  simulated  alone.  Dotted  line:  Drain  electrode  of  a  second 
operating  transistor  is  0.5  jam  apart  from  the  drain  of  the  simulated  transistor. 

B.  Electromagnetic-Physical  Modeling  of  Multifinger  Transistors 

It  is  clear  from  the  previous  results  that  EM-wave  propagation  and  electron-wave  interactions 
change  the  device  characteristics  at  high  frequency.  Accordingly,  different  structure  shapes  and 
configurations  need  to  be  employed  to  minimize  these  effects,  aiming  to  improve  the  device 
performance,  especially  at  high  operating  power  and  frequency.  A  possible  solution  is  to  use 
multiple  gate-fingers  of  shorter  lengths.  In  this  manner,  EM-wave  propagation  effects  along  the 
device  width  are  minimized.  Moreover,  attenuation  is  reduced  as  a  result  of  reducing  the  gate 
metallic  resistance.  However,  large  number  of  fingers  means  that  the  attenuation  and  EM-wave 
propagation  effects  are  increased  along  the  feeding  line.  Besides,  more  fingers  may  cause  more 
EM-waves  interference.  Thus,  EM-wave  synchronization  for  the  multiple  fingers  is  crucial  for 
maximum  power  and  minimum  interference.  Therefore,  the  number  of  fingers  and  distance 
between  gate-fingers  should  be  optimized  simultaneously.  It  is  noteworthy  to  say  that  EM-wave 
phase-velocity  mismatches  is  due  to  the  different  applied  voltages  to  the  electrodes  and  also  due 
to  the  unsymmetrical  shape  of  the  structure. 

Moreover,  for  the  case  of  the  four-finger  transistor,  the  shape  and  size  of  the  air-bridge 
connecting  the  different  fingers  affect  the  high-frequency  characteristics  of  the  transistor. 
Considering  Fig.  15.c,  it  should  be  noticed  that  new  capacitances  between  the  air-bridge  and 
transistor  electrodes  Cair  bridge  are  created.  This  would  definitely  change  the  EM-wave  phase 

velocities  and  as  a  result  change  the  device  behavior.  Thus  an  optimal  air-bridge  structure  and 
size  should  be  employed  as  well.  Furthermore,  the  air-bridge  should  not  be  fragile  in  order  not  to 
break  easily,  which  represents  an  extra  constraint  that  needs  to  be  included  in  our  optimization 
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problem.  The  feeding  line  shape  represents  also  a  parameter  that  needs  to  be  considered  for 
circuit-matching  issues. 

In  this  paper,  ad-hoc  optimization  is  performed  to  obtain  near-optimal  transistor  parameters 
based  on  the  above  criteria.  Table  IV  shows  the  new  parameters  for  the  optimized  multi  linger 
transistors,  and  Fig.  15  gives  a  generic  3D  view  of  the  simulated  multifinger  transistors. 

Output  voltages  for  the  simulated  multifinger  transistors  are  shown  in  Fig.  16.  Considering 
this  figure,  one  should  observe  that  the  voltage-gain  increases  when  using  four-finger  transistors. 
In  addition,  the  shape  of  the  output  signal  for  the  four-finger  transistor  case  appears  to  be  much 
better,  which  means  fewer  harmonics. 

Design  and  optimization  of  high-frequency  multifinger  transistors  need  a  tremendous 
research  work.  It  requires,  as  a  backbone,  a  very  efficient  numerical  simulator  that  includes  EM- 
wave  propagation  and  electron-wave  interactions.  This  simulator  should  be  accurate  and  most 
importantly  fast  in  order  to  be  suitable  for  optimizing  complex  microwave  structures.  It  is  our 
belief  that  this  paper  presents  this  type  of  simulator.  This  paper  also  presents,  for  the  first  time  in 
literature,  the  preliminary  numerical  results  of  electromagnetic  physical  simulation  of  multifinger 
transistors,  based  on  ad-hoc  optimization.  The  future  research  work  will  employ  rigorous 
optimization  techniques  to  obtain  the  optimal  multifinger  transistor  structure  based  on  the  model 
presented  in  this  paper.  Moreover,  measurements  will  be  carried  out  and  compared  to  the  results 
achieved  by  our  model. 


TABLE  IV 

Multifinger  Transistor  Optimized-Parameters  Used  In  The  Simulation 


Drain  and  source  contacts 
Gate-source  separation 
Gate-drain  separation 
Device  thickness 
Device  length 
Gate  length 
Device  Width 
Active  layer  thickness 
Active  layer  doping 
Schottky  barrier  height 
DC  gate-source  voltage 
DC  drain-source  voltage 
Operating  frequency 


0.5  pm 
0.5  pm 
0.4  pm 
0.4  pm 
2.1  pm 
0.2  pm 

1x450,  2x225,4x112.5  pm 
0.1  pm 
2xl017  cm"3 
0.8  V 
-0.2  V 
3.0  V 
60  GHz 
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Fig.  15.  Generic  3D  view  of  the  simulated  multifmger  transistors  (not  to  scale). 

(a)  Single-finger  transistor  (lx  450|am).  (b)  Two-finger  transistor  (2x  225|um). 
(c)  Four-finger  transistor  (4x1 12.5|um). 
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Fig.  16.  Output  voltage  for  the  simulated  multifmger  transistors  when 
EM-wave  propagation  and  electron- wave  interactions  are  considered. 

X.  Conclusions 

In  this  paper,  the  potential  of  high  power  and  frequency  multifinger  transistors  is 
demonstrated,  using  a  new  wavelet-based  full-wave  physical  simulation  approach.  The  proposed 
technique  solves  the  model  that  combines  the  transport  physics,  and  Maxwell’s  Equations  on 
nonuniform  self-adaptive  grids,  obtained  using  wavelets.  Moreover,  efficient  grid-updating 
criteria  for  the  active  and  passive  parts  of  the  problem  are  developed  and  verified.  A  reduction  of 
75%  in  CPU-time  is  achieved  compared  to  a  uniform  grid  case  with  an  error  of  2%  on  the  DC 
drain  current.  In  addition,  an  80%  CPU-time  reduction  is  obtained  for  FDTD  simulations  with 
approximately  0.1%  average  error  on  the  potential.  Tradeoffs  are  observed  between  the  threshold 
value,  CPU-time,  and  accuracy,  suggesting  an  optimal  value  for  the  threshold. 

The  preliminary  results  of  this  paper  show  that  at  very  high  frequency,  several  phenomena 
with  strong  impact  on  the  device  behavior  start  to  emerge,  such  as  phase  velocity  mismatches, 
electron-wave  interactions,  and  attenuation.  The  results  suggest  that  contemporary  microwave 
devices  should  be  optimized  to  minimize  these  effects  or  possibly  take  advantage  of  in  favor  of 
improved  device  characteristics.  The  results  also  recommend  multifinger  transistors  as  potential 
alternatives  to  conventional  transistors.  This  is  achieved  by  using  multiple-finger  gates  of  less 
width  instead  of  a  single-gate  device.  Furthermore,  this  paper  underlines  the  enhanced 
microwave  characteristics  of  multifinger  transistors  attributable  to  reducing  attenuation  and  EM- 
wave  propagation  effects  along  the  device  width.  The  future  research  work  will  involve 
employing  rigorous  optimization  techniques  to  obtain  the  optimal  multifinger  transistor  structure 
based  on  the  electromagnetic-physical  model  presented  in  this  paper.  Moreover,  measurements 
will  be  carried  out  and  compared  to  the  results  achieved  by  our  model. 
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Abstract —  This  paper  presents  a  new  approach  for  the  simulation  and  optimization 
of  microwave  devices,  using  a  genetic  algorithm  (GA).  The  proposed  technique  solves  the 
equations  that  describe  the  semiconductor  transport  physics  in  conjunction  with  Poisson’s 
equation,  employing  an  adaptive  real-coded  GA.  An  objective  function  is  formulated,  and 
most  of  the  GA  parameters  are  recommended  to  change  during  the  simulation.  In  addition, 
different  methods  for  describing  the  way  the  GA  parameters  change  are  developed  and 
studied.  The  effect  of  GA  parameters  including  the  mutation  value,  number  of  crossover 
points,  selection  criteria,  size  of  population,  and  probability  of  mutation  is  analyzed.  The 
technique  is  validated  by  simulating  a  submicrometer  field  effect  transistor  (FET),  and 
then  compared  to  successive  over  relaxation  (SOR),  showing  the  same  degree  of  accuracy 
along  with  a  moderate  speed  of  convergence.  The  purpose  of  this  paper  is  to  introduce  a 
new  vision  for  a  genetic  algorithm  capable  of  optimizing  real  value  functions  with  a 
considerably  large  number  of  variables.  This  work  also  represents  a  fundamental  step 
toward  applying  GA’s  to  Maxwell’s  equations  in  conjunction  with  the  hydrodynamic 
model  (HDM),  aiming  to  develop  an  optimized  and  unconditionally  stable  global-modeling 
simulator. 


I.  Introduction 

CONTEMPORARY  high  performance  electronics  are  based  on  technologies  such  as 
monolithic  microwave  integrated  circuits  (MMIC’s),  with  a  large  number  of  closely  packed 
passive  and  active  structures,  several  levels  of  transmission  lines  and  discontinuities,  all 
operating  at  high  speeds,  frequencies,  and  sometimes  over  very  broad  bandwidths.  It  is  thus 
anticipated  that  the  design  of  MMIC’s  should  involve  robust  design  tools  that  would  simulate  all 
the  circuit  elements  simultaneously.  The  possibility  of  achieving  this  type  of  modeling  is 
addressed  by  global  circuit  modeling  that  has  been  demonstrated  in  [l]-[4]. 
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Global  modeling  is  a  tremendous  task  that  involves  advanced  numerical  techniques,  and 
different  algorithms  with  tight  stability  constraints  [3].  Thus,  there  is  an  urgent  need  to  develop 
and  present  new  simulation  approaches  that  could  relax  or  even  eliminate  these  constraints.  On 
the  other  hand,  genetic  algorithms  are  numerical  optimization  algorithms  inspired  by  both  natural 
selection  and  natural  genetics,  which  are  unconditionally  stable.  The  method  is  a  general  one, 
capable  of  being  applied  to  an  extremely  wide  range  of  problems.  GA’s  have  proven  themselves 
for  optimizing  many  large  and  complex  problems  in  our  field  [5]-[l  1  ]. 

Some  optimization  problems  have  multiple  local  minima.  Where  methods  based  on  steepest 
descent  would  fall  in  one  of  these  local  minima,  resulting  in  a  different  solution.  GA’s  are 
random  algorithms  and  researchers  have  found  their  generality  and  that  they  are  unconditionally 
stable.  GA’s  are  thus  suitable  to  find  the  global  solution  for  problems  having  multiple  minima 
[12].  Furthermore,  in  many  problems  requiring  solving  systems  of  linear  equations  Ax  =  b ,  the 
matrix  A  has  a  large  Condition  Number.  For  these  problems,  standard  methods  will  not  be  able 
to  get  the  correct  solution.  For  instance,  solving  Poisson’s  equation  on  a  nonuniform  grid.  In 
this  case,  genetic-based  algorithms  would  outperfonn  standard  methods.  This  is  because  standard 
methods  work  fine  only  for  well-posed  problems  (problems  with  A  having  a  small  condition 
number).  On  the  other  hand,  a  genetic -based  algorithm  converges  independent  of  the  condition 
number  of  A . 

From  the  above,  it  is  motivating  to  make  an  effort  to  apply  GA’s  to  MaxwelFs  equations  or 
the  hydrodynamic  model  (HDM),  aiming  to  develop  an  optimized  and  unconditionally  stable 
algorithm.  It  is  noteworthy  to  say  that  the  main  purpose  of  this  paper  is  to  lay  the  foundation  of  a 
genetic  algorithm  capable  of  optimizing  real  value  problems,  with  a  considerably  large  number 
of  unknowns. 

In  this  stage  of  the  work,  we  will  demonstrate  that  genetic  algorithms  can  be  applied  to  the 
hydrodynamic  model  (HDM)  in  conjunction  with  Poisson’s  equation  to  accurately  model 
submicrometer  gate  devices,  with  less  stability  constraints.  Ultimately,  a  hydrodynamic  model 
should  be  implemented  with  equations  that  would  have  numerical  stability  restrictions  such  as 
MaxwelFs  equations  rather  than  Poisson’s  equation  in  order  to  obtain  a  self-consistent 
simulation  of  electromagnetic-wave  propagation  effects,  employing  an  optimized  and 
unconditionally  stable  algorithm. 

This  paper  is  organized  as  follows.  Section  II  describes  the  problem  under  consideration. 
While,  section  III  presents  in  details  the  implementation  of  the  proposed  algorithm.  Results 
along  with  illustrative  graphs  are  provided  in  section  IV.  Finally,  conclusions  are  presented  in 
section  V. 


II.  Problem  Description 
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The  transistor  model  used  in  this  work  is  a  two  dimensional  full-hydrodynamic  large-signal 
physical  model.  The  active  device  model  is  based  on  the  moments  of  the  Boltzmann’s  Transport 
equation  obtained  by  integrating  over  the  momentum  space.  The  integration  results  in  a  strongly 
coupled  highly  nonlinear  set  of  partial  differential  equations,  called  the  conservation  equations. 
These  equations  provide  a  time-dependent  self-consistent  solution  for  carrier  density,  carrier 
energy,  and  carrier  momentum,  which  are  given  as  follows. 

■  current  continuity 

■  energy  conservation 


—  +  V.(no)  =  0. 
dt 


(1) 


d(ne) 

dt 


+  qnv.E  +  V.(«u( s  +  KBT)  )  =  - 


»(s-s0) 

TsOO 


(2) 


■  x-momentum  conservation 


d(^x)  +  WEX  +  (npx  or  +  nKBT  ) 

dt  ox 


<PxzlA 

^»oo 


(3) 


In  the  above  equations,  n  is  the  electron  concentration,  o  is  the  electron  velocity,  E  is  the 
electric  field,  s  is  the  electron  energy,  s0  is  the  equilibrium  thennal  energy,  and  p  is  the  electron 
momentum.  The  energy  and  momentum  relaxation  times  are  given  by  r£  and  rm ;  respectively. 
Similar  expression  can  be  obtained  for  the  y-direction  momentum.  The  three  conservation 
equations  are  solved  in  conjunction  with  Poisson’s  equation: 


V2<P  =  ^(Nd-n)  (4) 

where  (j)  is  the  electrostatic  potential,  q  is  the  electron  charge,  e  is  the  dielectric  constant,  N d  is 
the  doping  concentration,  and  n  is  the  carrier  concentration  at  any  given  time.  The  total  current 
density  distribution  J  inside  the  active  device  at  any  time  t  is  given  as: 


J(t)  =  -qno(t). 


The  low  field  mobility  is  given  by  the  empirical  relation  [13]: 


Mo=- 


8000 


1  + 


K-io-'T 


cm"  /  V.s. 


(5) 

(6) 
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On  the  other  hand,  the  mobility  for  large-signal  simulations  is  calculated  as  ju  =  od  /  Ess , 
where  vd  is  estimated  using  Eq.  7. 

'  E„ 


Mo  Ess  +  VS! 


0,1  = 


4500 


1  + 


4500 


(V) 


In  the  above  equation,  vd  is  the  electron  drift-velocity,  ju0  is  low  field  mobility  given  by  Eq. 
6,  oss  and  Ess  are  the  steady-state  electron  velocity  and  electric  field,  respectively.  It  is 
significant  to  note  that  both  uss  and  Ess  are  functions  of  energy,  and  they  get  updated  each  time  a 
new  energy  distribution  is  estimated  using  the  hydrodynamic  model. 

The  above  model  accurately  describes  all  the  non-stationary  transport  effects  by 
incorporating  energy  dependence  into  all  the  transport  parameters  such  as  the  effective  mass  and 
relaxation  times.  Fig.  1  shows  the  cross-section  of  the  simulated  structure. 
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Fig.  1 .  Cross-section  of  the  simulated  transistor. 


III.  Optimization  Using  Genetic  Algorithms 

In  this  section,  we  apply  a  genetic  algorithm  for  the  solution  of  the  boundary  value 
problem  for  the  distribution  of  potential  across  an  FET. 

A  generic  flowchart  of  the  algorithm  is  shown  in  Fig.  2.  The  first  step  is  to  read  the 
matrix  A  and  vector  b  of  the  system  of  linear  equations,  Ax  =  b ,  which  are  derived  from 
Poisson’s  equation  [14].  A  population  of  random  solutions  (chromosomes)  representing  the 
vector  solution  x  is  then  initialized.  Next,  the  objective  function  (fitness)  is  estimated  for  all  the 
chromosomes  that  have  been  randomly  generated.  Based  on  the  fitness,  two  parents  are 
generated  either  by  roulette  wheel  selection  or  by  tournament  selection  methods.  Mutation  and 
crossover  are  then  performed  on  the  selected  parents  to  generate  two  children.  Replacement  is 
conducted  by  comparing  the  fitness  of  children  with  their  parents,  and  the  worst  two 
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chromosomes  are  removed  from  the  population.  The  best  chromosome  is  identified  based  on 
fitness,  and  finally  a  check  is  carried  out  against  a  certain  stopping  criteria  to  either  stop  the 
simulator,  or  to  perform  another  iteration.  The  details  of  implementing  the  proposed  algorithm 
are  as  follows. 


Fig. 2.  Generic  flowchart  of  the  genetic  algorithm 


•  Step  One:  Initialization 
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Initialization  is  done  by  randomly  generating  M  chromosomes  representing  the  GA 
population.  Real  encoding  is  adopted  for  our  problem.  Each  chromosome  contains  n  genes, 
which  corresponds  to  the  variables  in  the  vector  solution  x .  The  generated  random  numbers 
have  a  range  associated  with  the  applied  DC  voltages  to  the  device  electrodes. 


•  Step  Two:  Evaluate  Fitness 

Each  chromosome  is  evaluated  based  on  an  objective  function.  The  objective  function  is 
developed  in  a  way  that  it  accurately  determines  how  close  the  randomly  generated  solutions  are 
to  the  optimal  solution. 

Dealing  with  2D  Poisson’s  equation  means  that  the  matrix  A  has  one  to  five  elements  in 
each  row  (sparse  matrix)  [14].  Based  on  this,  it  is  found  that  the  vector  solution  x  with  size 
n  =  Nx  N  v  can  be  fully  constructed  by  knowing  only  the  2NX  elements  close  to  the  electrodes 

following  the  flowchart  given  in  Fig.  3.  Fig.  4  shows  a  five-by-five-grid  to  illustrate  the 
implementation  of  the  algorithm  shown  in  Fig.  3.  Each  square  in  Fig.  4  represents  a  grid-point.  In 
Fig.  4,  the  elements  inside  the  gray  squares  are  used  to  estimate  the  next  element,  which  is  inside 

the  darker  square.  For  instance,  x23  is  estimated  as  a  linear  combination  of  xl9,  x24 ,  and  x25 , 
following  the  algorithm  in  Fig. 3.  This  process  is  repeated  until  all  elements  in  our  domain  are 
estimated,  as  shown  in  Fig.4.  It  is  worth  mentioning  that  the  proposed  numbering  sequence  is 
crucial  for  the  correct  estimation.  On  the  other  hand,  the  2 Nx  elements  needed  to  implement  the 

algorithm  in  Fig.  3  can  be  randomly  generated  with  a  minimal  error  based  on  the  following. 
Since  boundary  conditions  must  be  satisfied,  then  the  values  of  the  potential  are  known  precisely 
at  the  boundaries.  Moreover,  an  estimate  for  the  value  of  the  potential  near  boundaries  can  be 
randomly  generated  bearing  in  mind  that  they  would  have  very  close  values  to  the  potential  at  the 
electrodes.  The  randomly  generated  solution  g  along  with  the  second  nonn  of  Ax-b  are 
included  in  the  objective  function  given  by  Eq.  8,  which  needs  to  be  maximized  for  a  minimum 
value  of  the  dominator.  In  Eq.  8,  g  is  the  randomly  generated  solution,  while  w  is  a  scalar 
representing  the  weight  of  the  nonn  criterion.  The  value  of  w  is  chosen  to  be  10%,  by  trial  and 
error. 


Where, 


Fitness  - 


I  Zc, 


1  +  w  '  ,  =  l - +  D 

n 


uj"'  ]  "i 


(8) 
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Fig. 3.  Flowchart  of  the  randomly  generated  solution  G. 
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•  Step  Three:  Regeneration 

Two  methods  are  used  for  parent  regeneration,  namely,  Roulette  Wheel  and  Tournament 
selections.  The  details  of  implementing  each  method  are  described  below. 

A.  Roulette  Wheel  Selection 

Parents  in  the  Roulette  Wheel  method  are  chosen  randomly  according  to  their  fitness.  As 
the  name  implies,  the  method  imitates  the  Roulette  Wheel  game,  where  the  thrown  dice  would 
most  probably  end  being  in  the  slot  with  the  largest  area.  Following  this,  one  can  conclude  that 
the  chromosome  with  the  largest  fitness  value  is  most  likely  to  be  chosen  because  it  has  the 
largest  slot  size  or  equivalently  the  largest  fitness. 

B.  Tournament  Selection 

In  this  method,  two  groups  from  the  population  are  randomly  selected,  sub-populations.  It  is 
worth  mentioning  that  the  population  size  is  chosen  randomly  as  well.  The  best  chromosome 
from  each  of  the  randomly  generated  sub-populations  is  chosen  to  represent  a  parent. 

Step  Four:  Crossover 

Now,  two  parents  have  been  selected  for  their  genes  to  be  crossed  over  and  mutated. 
Crossover  is  conducted  by  first  randomly  selecting  a  crossover  point  within  the  chromosome. 
Two  children  are  then  conceived  by  mixing  the  genes  of  the  two  parents  at  the  crossover  point. 
At  this  moment,  two  different  parameters  can  be  analyzed.  The  first  parameter  is  the  number  of 
crossover  points,  i.e.,  more  than  one  crossover  point  can  be  achieved.  The  other  parameter  is  the 
number  of  genes  involved  in  each  crossover  point,  and  will  be  denoted  by  the  crossover  width. 
The  effects  of  both  parameters  are  studied  and  included  in  the  results  section. 

•  Step  Five:  Mutation 

Mutation  is  carried  out  by  randomly  changing  one  or  more  genes  (variables)  of  the 
created  offspring.  We  then  have  two  mutation  parameters  to  study  their  effect.  The  first  one  is 
the  number  of  mutated  genes  or  variables  within  the  chromosome.  The  other  parameter  deals 
with  the  value  of  mutation. 

Step  Six:  Replacement 

Replacement  is  performed  by  comparing  the  fitness  of  the  parents  with  their  offspring.  The 
best  two  chromosomes  out  of  the  four  are  included  in  the  population  for  the  next  iteration. 

•  Step  Seven:  Ending  the  Algorithm 

The  best  chromosome  is  identified  at  each  iteration  and  error  is  calculated  as  the  second 
norm  of  Ax-b  .  This  error  is  checked  against  a  predefined  value,  if  satisfied,  the  simulator  stops 
and  prints  the  final  results. 
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IV.  Results  And  Discussions 


In  this  section,  the  effect  of  different  GA  parameters  on  the  algorithm  behavior  is 
investigated.  The  GA  parameters  used  in  the  simulation  are  summarized  in  Table  I. 


Table  I 

GA  Parameters  Used  In  The  Simulation 


MU 

Mutation  factor  or  value 

NC 

Number  of  crossover  points 

CW 

Crossover  width 

NM 

Number  of  mutated  variables 

SC 

Selection  criterion 

NPOP 

Size  of  population 

PM 

Probability  of  mutation 

The  default  values  are  0.1,  1,  1,  1,  1,  100,  and  1  for  MU,  NC,  CW,  NM,  SC,  NPOP,  and 
PM,  respectively. 

It  is  significant  to  note  here  that  the  developed  algorithm  is  implemented  as  a  subroutine 
to  solve  the  system  of  linear  equations,  Ax  =  b .  Poisson’s  equation  is  then  coupled  to  the  HDM 
equations  as  a  subroutine.  The  coupling  is  carried  out  as  follows.  First,  the  hydrodynamic 
equations  are  solved  to  get  the  updated  value  for  carrier  density  n  .  The  updated  carrier  density  is 
then  plugged  into  Poisson’s  equation,  resulting  in  a  new  system  of  linear  equations.  The  new 
system  of  linear  equations  is  passed  to  the  genetic-based  Poisson  solver  to  solve  for  x ,  i.e.,  the 
updated  value  of  the  potential.  The  potential  is  differentiated  to  get  the  electric  filed.  Finally,  the 
updated  value  of  the  electric  field  is  plugged  into  the  HDM  to  estimate  the  updated  value  of 
carrier  density.  This  process  is  repeated  until  the  stopping  criterion  is  satisfied. 

Fig.  5  shows  the  distance  from  the  optimal  solution  versus  number  of  generations  for 
different  values  of  MU.  The  mutation  value  of  any  gene  (variable)  is  proposed  as  follows. 

ChildneJi)  =  Childold(i)  +  MU  ■  RND  (9) 

Where  MU  is  the  mutation  factor  and  RND  is  a  random  number  between  -1  and  1.  Fig.  5  shows 
that  the  best  result  is  obtained  when  the  mutation  value  is  dependent  on  the  fitness.  The  reason  is 
that  as  the  value  of  the  fitness  increases,  which  indicates  being  very  close  to  the  optimal  solution, 
the  mutation  factor  decreases.  In  this  manner,  the  mutation  value  is  changed  in  the  correct  way 
for  a  faster  convergence.  Moreover,  introducing  a  random  feature  along  with  the  dependence  of 
MU  on  the  fitness  does  not  enhance  the  convergence.  A  general  conclusion  is  that  smaller  values 
of  MU  are  observed  to  have  better  convergence  curves. 
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Figures  6  and  7  illustrate  the  effect  of  the  number  of  crossover  points  and  crossover 
width,  respectively.  From  Fig.  6,  it  can  be  concluded  that  larger  number  of  crossover  points 
is  the  right  choice  for  better  accuracy  along  with  higher  speed  of  convergence.  Moreover, 
choosing  the  number  of  crossover  points  to  change  randomly  within  simulation  does  not 
improve  the  algorithm.  On  the  other  hand,  considering  Fig.  7,  it  is  apparent  that  rate  of 
convergence  of  the  genetic  algorithm  is  independent  of  the  crossover  width  value. 


Fig.  5.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  mutation  values. 

Fig.  8  shows  the  effect  of  the  number  of  mutated  elements  on  the  convergence  and 
accuracy  of  the  proposed  algorithm.  It  can  be  observed  that  as  the  number  of  the  mutated 
variables  decreases,  the  convergence  and  accuracy  of  the  algorithm  are  improved.  The  best 
curve  is  obtained  for  a  number  of  mutated  variables  equals  to  one  percent.  This  complies 
with  nature,  since  biological  mutation  hits  only  a  very  small  number  of  genes.  Moreover, 
changing  the  number  of  the  mutated  variables  randomly  throughout  the  simulation  introduces 
a  reasonable  improvement  over  the  one  percent  mutation  case. 
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Fig.  6.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  number  of  crossover  points. 
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Fig.  7.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  values  of  crossover  widths. 


Fig.  9  demonstrates  how  the  choice  of  parent  selection  method  affects  the  algorithm 
performance.  It  can  be  pointed  out  that  roulette  wheel  selection  has  the  best  performance. 
Furthermore,  employing  a  hybrid  technique  does  not  improve  the  algorithm.  For  instance,  error 

_2 

reaches  10  in  almost  500  generations  when  roulette  wheel  selection  is  employed,  whereas 
1200  generations  are  needed  for  tournament  selection  to  reach  the  same  value  of  error.  It  is 
worth  mentioning  here  that  roulette  wheel  selection  inherently  uses  some  sort  of  elitism. 
Employing  elitism  may  or  may  not  be  useful  depending  on  the  problem  under  consideration.  The 
main  reason  for  roulette  wheel  selection  producing  better  results  is  the  choice  of  the  objective 
function  given  by  Eq.  8.  It  is  important  to  mention  that  tournament  selection  is  known  to  produce 
better  results  over  the  roulette  wheel  method.  However,  this  is  not  general,  and  the  numerical 
example  provided  in  this  paper  emphasizes  that. 


Fig.  8.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  numbers  of  mutated  elements. 

Fig.  10  shows  the  effect  of  the  probability  of  mutation.  It  should  be  noticed  that 
increasing  the  probability  of  mutation  has  a  positive  effect  on  the  perfonnance  of  the  algorithm, 
and  the  mutation  probability  is  not  crucial  to  the  algorithm  as  long  as  it  is  relatively  high.  Fig.  1 1 
shows  the  effect  of  the  population  size  on  the  algorithm  convergence.  This  figure  emphasizes 
that  the  population  size  is  not  a  critical  parameter.  The  reason  is  that  the  proposed  objective 
function  given  by  Eq.  8,  inherently  allows  elitism.  This  makes  the  proposed  genetic  algorithm 
independent  of  population  size. 
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Fig.  9.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  selection  criteria. 


number  of  generations 

Fig.  10.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  probabilities  of  mutation. 
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Fig.  11.  Distance  from  the  optimal  solution  versus  number  of  generations 
for  different  population  sizes. 

A.  DC  Simulation  Results 

To  demonstrate  the  potential  of  the  proposed  approach,  it  is  applied  to  an  idealized 
MESFET  structure,  which  is  discretized  by  a  mesh  of  32  Ac  by  32  Av  with  At  =  0.001  ps . 
Forward  Euler  is  adopted  as  an  explicit  finite  difference  method.  In  addition,  upwinding  is 
employed  to  have  a  stable  scheme.  The  space  step  sizes  are  adjusted  to  satisfy  Debye  length, 
while  the  time  step  value  At  is  chosen  to  satisfy  the  Courant-Friedrichs-Levy  (CFL)  condition. 
While,  Poisson’s  equation  is  solved  using  the  proposed  algorithm. 

Fig.  12. a  shows  the  potential  distribution  obtained  using  the  proposed  algorithm.  This 
graph  demonstrates  that  boundary  conditions  are  satisfied  at  the  electrodes.  For  instance,  the 
value  of  the  potential  at  the  gate  equals  to  -1.3  volts,  which  is  the  applied  DC  voltage  minus  the 
Schottky  barrier  height.  While,  Fig.  12. b  shows  the  carrier  density  distribution.  It  is  significant 
to  indicate  that  the  proposed  algorithm  gives  precisely  the  same  results  obtained  when  SOR 
algorithm  is  employed.  The  comparison  results  between  the  algorithms  are  not  provided  because 
their  results  coincide  exactly  on  each  other.  It  is  noteworthy  to  say  that  the  purpose  of  this 
section  is  to  show  that  genetic  algorithms  can  be  applied  to  solve  real  value  problems  having  a 
large  number  of  unknowns,  with  a  very  high  degree  of  accuracy.  The  speed  of  convergence  is  not 
an  objective  at  this  stage  of  the  work.  Ultimately,  the  developed  genetic  algorithm  needs  to  be 
applied  to  equations  that  have  stability  constraints  in  order  to  have  unconditionally  stable 
algorithm,  or  to  solve  problems  that  traditional  optimization  techniques  cannot  solve  (problems 
with  multiple  local  minima). 
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Fig.  12.  Sample  DC  Results  obtained  using  the  proposed  algorithm, 
(a)  Potential  distribution,  (b)  Carrier  density  distribution. 
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B.  AC  Simulation  Results 

The  AC  excitation  applied  to  the  gate  electrode  is  given  as: 

Vgs(t)  =  Vgs0  +  A  vgssin{cot)  (10) 

where  Vgs  is  the  DC  bias  applied  to  the  gate  electrode,  A  ogs  is  the  peak  value  of  the  AC  signal 
(0.1  volts),  and  a>  is  the  frequency  of  the  applied  signal  in  rad./ sec.  The  frequency  used  in  the 
simulation  is  60  GHz. 

First,  the  DC  distribution  is  obtained  by  solving  Poisson’s  equation,  using  the  proposed 
algorithm  in  conjunction  with  the  three  hydrodynamic  conservation  equations.  Then,  a  new  value 
of  Vgs  is  calculated  using  Eq.  10.  The  new  value  of  Vgs  is  used  to  update  Poisson’s  equation  to 
get  the  new  potential  distribution.  The  electric  field  is  then  estimated  and  used  to  update  the 
variables  in  the  conservation  equations.  This  process  is  repeated  every  At  until  t  =  tmax .  The 
current  density  is  obtained  using  Eq.  5.  The  current  density  calculated  on  the  plan  located 
midway  between  the  drain  and  gate  is  integrated  to  obtain  the  total  current.  The  output  voltage  is 
estimated  by  multiplying  the  total  current  by  the  resistance  that  defines  the  DC  operating  point 
(Q  point)  of  the  transistor.  Fig.  13  shows  the  AC  gate  and  drain  voltages.  A  gain  of  11  dB  is 
achieved.  Moreover,  it  is  observed  that  there  is  an  output  delay  of  about  1  ps  that  represents  the 
time  required  for  the  transistor  to  respond  to  the  input  signal. 


Fig.  13.  AC  gate  and  drain  voltages  obtained  using  the  proposed  algorithm. 
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V.  Conclusions 

In  this  paper,  a  new  technique  is  developed  for  solving  the  hydrodynamic  model  (HDM) 
in  conjunction  with  Poisson’s  equation,  using  an  adaptive  real-coded  genetic  algorithm.  Several 
GA  design  parameters  have  been  studied  to  illustrate  their  effects  on  the  algorithm  convergence. 
The  novelty  of  the  proposed  technique  comes  from  the  genetic  algorithm  itself.  This  has  been 
achieved  by  developing  a  very  efficient  objective  function,  along  with  introducing  completely 
new  concepts  such  as  fitness-dependent  GA  parameters.  Moreover,  the  problem  this  paper 
presents  is  a  new  application  for  genetic  algorithms.  In  addition,  the  proposed  genetic  algorithm 
will  outperform  standard  methods  for  several  types  of  problems.  For  instance,  finding  the  global 
solution  of  optimization  problems  having  multiple  local  minima,  and  problems  where  matrices 
have  large  Condition  Numbers.  This  paper  also  represents  a  fundamental  step  toward  applying 
GA’s  to  Maxwell’s  equations  in  conjunction  with  the  HDM,  aiming  to  develop  an  optimized  and 
unconditionally  stable  global-modeling  simulator. 
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